In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

import igrins_mod as ig # Custom module file for igrins shenanigans

import glob
import os
import warnings

import astropy.units as u
from astroquery.nist import Nist # atomic lines
# from astroquery.linelists.cdms import CDMS # molecular lines?

from lmfit import Model, Parameters, CompositeModel
from lmfit.models import GaussianModel

from scipy.stats import chisquare
from scipy.integrate import trapz, simpson
from scipy.optimize import curve_fit

from astropy.io import fits

# Plotting Parameters
plt.rcParams['font.size'] = 20
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 18
plt.rcParams['ytick.labelsize'] =18

plt.rcParams['legend.fontsize'] = 16
plt.rcParams['figure.titlesize'] = 20

plt.rcParams['axes.labelweight']='bold'
plt.rcParams['axes.linewidth'] = 3

plt.rcParams['xtick.major.size'] = 10
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['xtick.minor.size'] = 5

plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams['ytick.minor.size'] = 5

plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'

%matplotlib inline

Things that should only have to be defined once, here at the top¶


In [ ]:
# Size of 1 spectral resolution element
# IGRINS Spectral Resolution
spec_res = 0.00001

c = 299792458  # speed of light m/s

# Use Normalized (single) Gaussian Distribution
def gaussian_func(x, ampl, center, std):
    return ((ampl) / (std * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - center) / std) ** 2)) + 1

# Reduced and order-merged data filepath 
# Laptop Path
data_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-Spectra\\IGRINS_Merged"

# File path for figures to live in
fig_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-SpectraIGRINS_figs\\standards_spectra"

# Create the folder if it doesn't exist
if not os.path.exists(fig_path):
    os.makedirs(fig_path)

# Nicole's merged K-band spectra of some Taurus Standards
merged_standard_files = glob.glob(data_path + "/merged_standards/m*.fits")
standard_table = pd.read_csv('./standard_table_v2.txt', index_col=0)  # csv of standards with file and Spectral Type, c/v TBA

# early_k = ('1', '2', '3', '4', '5')
# Only want to look at the K types for the moment
# standard_table_k = standard_table[standard_table['Spectral_Type'].str.startswith(('K1', 'K2', 'K5'))]

standard_shift = standard_table['Wavelength Shift'].values
standard_list = standard_table['Source'].values

hops_table = pd.read_csv('./hops_table.txt')
hops_list = hops_table['Source']

# Determine the maximum length of flux arrays for the standards
max_flux_length = max(len(fits.getdata(file)[1]) for file in standard_list)
max_wavelen_length = max(len(fits.getdata(file)[0]) for file in standard_list)
max_snr_length = max(len(fits.getdata(file)[2]) for file in standard_list)

hops_flux_length = max(len(fits.getdata(file)[1]) for file in hops_table['Source'])
hops_wavelen_length = max(len(fits.getdata(file)[0]) for file in hops_table['Source'])
hops_snr_length = max(len(fits.getdata(file)[2]) for file in hops_table['Source'])

# Initialize stacks with NaN values
wavelen_stack = np.full((max_wavelen_length, len(standard_list)), np.nan)
raw_flux_stack = np.full((max_flux_length, len(standard_list)), np.nan)
snr_stack = np.full((max_snr_length, len(standard_list)), np.nan)
raw_flux_err_stack = np.full((max_snr_length, len(standard_list)), np.nan)

# hops_wavelen_stack = np.full((max_wavelen_length, len(hops_table)), np.nan)
# hops_raw_flux_stack = np.full((max_flux_length, len(hops_table)), np.nan)
# hops_snr_stack = np.full((max_snr_length, len(hops_table)), np.nan)
# hops_raw_flux_err_stack = np.full((max_snr_length, len(hops_table)), np.nan)


# Fill stacks with data
for i, file in enumerate(standard_list):
    # Get wavelength, flux, snr per resolution element data
    data = fits.getdata(file)
    wavelen, flux, snr = data[0], data[1], data[2]
    
    # Clean data a bit
    snr_min = 50  # Minimum SNR
    snr_max = 1e4  # Maximum SNR
    snr_cut = (snr > snr_min)  # & (snr < snr_max)  # Bitwise SNR masking

    flux_min = 10000  # Minimum flux
    flux_cut = flux > flux_min  # Bitwise flux masking

    wavelen_min = 2.0
    wavelen_max = 2.4
    wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)

    # Apply masks and remove NaNs and infs
    mask = flux_cut & wavelen_cut & snr_cut
    
    wavelen = wavelen[mask]
    flux = flux[mask]
    snr = snr[mask]

    # Remove NaNs and infs from wavelen, flux, and snr arrays
    valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
    
    wavelen = wavelen[valid_indices] - standard_shift[i]
    flux = flux[valid_indices]
    snr = snr[valid_indices]

    # Check for NaNs in the final arrays
    if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
        print(f"NaNs found in data for file {file} after cleaning")

    wavelen_stack[:len(wavelen), i] = wavelen  # Wavelength arrays for each standard
    raw_flux_stack[:len(flux), i] = flux
    snr_stack[:len(snr), i] = snr
    raw_flux_err_stack[:len(flux), i] = flux / snr

# for i, file in enumerate(hops_list):
#     # Get wavelength, flux, snr per resolution element data
#     data = fits.getdata(file)
#     wavelen, flux, snr = data[0], data[1], data[2]
    
#     # Clean data a bit
#     snr_min = 50  # Minimum SNR
#     snr_max = 1e4  # Maximum SNR
#     snr_cut = (snr > snr_min)  # & (snr < snr_max)  # Bitwise SNR masking

#     flux_min = 10000  # Minimum flux
#     flux_cut = flux > flux_min  # Bitwise flux masking

#     wavelen_min = 2.0
#     wavelen_max = 2.4
#     wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)

#     # Apply masks and remove NaNs and infs
#     mask = flux_cut & wavelen_cut & snr_cut
    
#     wavelen = wavelen[mask]
#     flux = flux[mask]
#     snr = snr[mask]

#     # Remove NaNs and infs from wavelen, flux, and snr arrays
#     valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
    
#     wavelen = wavelen[valid_indices] - standard_shift[i]
#     flux = flux[valid_indices]
#     snr = snr[valid_indices]

#     # Check for NaNs in the final arrays
#     if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
#         print(f"NaNs found in data for file {file} after cleaning")

#     hops_wavelen_stack[:len(wavelen), i] = wavelen  # Wavelength arrays for each standard
#     hops_raw_flux_stack[:len(flux), i] = flux
#     hops_snr_stack[:len(snr), i] = snr
#     hops_raw_flux_err_stack[:len(flux), i] = flux / snr
In [ ]:
# Directly query NIST to find line features in K-band
with warnings.catch_warnings():  # Ignore warnings
    warnings.simplefilter('ignore') 
    lines_table = Nist.query(2.08*u.um,2.35*u.um,
                    linename = 'Na I, Sc I, Si I, Fe I, Fe II, Al I, Mg I, Ca I, H I, Ti I',
                    energy_level_unit='eV',output_order='wavelength')

igrins_wav_cut = (lines_table['Observed'] > 2.08) & (lines_table['Observed'] < 2.35)
lines_table = lines_table[igrins_wav_cut]
# lines_table = pd.read_csv('lines_table.txt')

# Make masks for the table of all the lines just in case I want to peek at certain transitions/wavelengths
na1_mask = lines_table['Spectrum'] == 'Na I'
sc1_mask = lines_table['Spectrum'] == 'Sc I'
si1_mask = lines_table['Spectrum'] == 'Si I'
fe1_mask = lines_table['Spectrum'] == 'Fe I'
fe2_mask = lines_table['Spectrum'] == 'Fe II'
al1_mask = lines_table['Spectrum'] == 'Al I'
mg1_mask = lines_table['Spectrum'] == 'Mg I'
ca1_mask = lines_table['Spectrum'] == 'Ca I'
h1_mask  = lines_table['Spectrum'] == 'H I'
ti1_mask = lines_table['Spectrum'] == 'Ti I'

# Just add all the masks to a list for the sake of my plotting a few cells down
mask_list = [na1_mask,sc1_mask,si1_mask,fe1_mask,al1_mask,mg1_mask,ca1_mask,h1_mask,ti1_mask]
color_list = ['purple', 'orange', 'green', 'blue', 'brown', 'crimson', 'olive', 'cyan', 'darkgreen']

Veiling (?) Lines¶

Fe I ~ 2.084 (lines_table[fe1_mask][6]), 2.275 (lines_table[fe1_mask][109]), 2.284 (lines_table[fe1_mask][115])

Si I ~ 2.092 (lines_table[si1_mask][0])

Mg I ~ 2.106 (lines_table[mg1_mask][0]), 2.282 lines_table[mg1_mask][11]

Al I ~ 2.11 (lines_table[al1_mask][1])

Ca I ~ 2.283

Ti I ~ 2.29 (lines_table[ti1_mask][69])

CO (2-0) > 2.294

In [ ]:
lines_table[fe1_mask][6]
Out[ ]:
Row index=6
SpectrumObservedRitzTransitionRel.AkifikAcc.Ei EkLower levelUpper levelTypeTPLine
str5float64str18float64str17float64float64str3str40str45str59str5str17str20
Fe I2.084652.08464934796.9681200------6.01718438 - 6.611932863d7.(4F).5s | e 3F | 33d7.(4F).5p | 3F* | 4----L11631

Fe I 2.08465 $\mu m$¶

In [ ]:
line_name = lines_table[fe1_mask][0]['Spectrum']
line_center = lines_table[fe1_mask][6]['Observed']

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)

    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(rf"{line_name} {line_center}$\mu m$ {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
amplitude1_stack = []
center1_stack = []
sigma_stack = []

pcov_stack = []

best_model_stack = []

flux_constant = np.linspace(0,-1,len(standard_table))

# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)


# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev = 50000)

    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma_stack.append(popt[2])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# for i in range(len(standard_table)):
    # fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
    # area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
    # ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    
fe1_2_0846_fit = []

for i in range(len(standard_table)):
    fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
    
    area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)

    # print(f'{area_int} {area_v2}')

ew1_stack = [] # empty list to load in equivalent widths

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))

    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
    

    # plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
    # Fe I
    # plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
    # plt.ylim(0)
    plt.ylabel('Flux + Constant')
    plt.xlabel(r'Wavelength [$\mu$m]')

    plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.legend(loc='best')
    plt.show()
In [ ]:
# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][109]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
amplitude1_stack = []
center1_stack = []
sigma_stack = []

pcov_stack = []

best_model_stack = []

flux_constant = np.linspace(0,-1,len(standard_list))

# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)


# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)

    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma_stack.append(popt[2])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# for i in range(len(standard_table)):
    # fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
    # area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
    # ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    
# fe1_2_0846_fit = []

# for i in range(len(standard_table)):
#     fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
    
#     area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
#     area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)

#     # print(f'{area_int} {area_v2}')

# ew1_stack = [] # empty list to load in equivalent widths

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))

    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
    

    # plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
    # Fe I
    # plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
    # plt.ylim(0)
    plt.ylabel('Flux + Constant')
    plt.xlabel(r'Wavelength [$\mu$m]')

    plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.legend(loc='best')
    plt.show()
c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

Si 1¶

In [ ]:
# Define the region for fitting
line_name = lines_table[si1_mask][0]['Spectrum'] # Species
line_center = lines_table[si1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {line_center} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Mg I¶

In [ ]:
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Al I¶

In [ ]:
# Define the region for fitting
line_name = lines_table[al1_mask][0]['Spectrum'] # Species
line_center = lines_table[al1_mask][1]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
 

Ti I¶

In [ ]:
# Define the region for fitting
line_name = lines_table[ti1_mask][0]['Spectrum'] # Species
line_center = lines_table[ti1_mask][69]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
 

Temperature Sensitive Lines¶

Sc I (2.2058 $\mu m$ and 2.2071 $\mu m$): index 1 and 2

Si I (2.2068 $\mu m$)

Mg/Al (~2.11 $\mu m$)

Fe I (2.2205–2.2346 $\mu m$) Over Ti region: index 10 and 11

Ca I (2.2614, 2.2631, 2.2657 $\mu m$): Only three lines

Na Interval ~2.21¶

Na I Doublet (2.2062, 2.2089) & Sc I (2.2058)

In [ ]:
# Define the region for fitting
line_name = lines_table[na1_mask][0]['Spectrum'] # Species
line_center = lines_table[na1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,30),(-140,30),(150,40),(350,30)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
# Initialize storage lists
params_stack = [] # list of best of fit parameters
params_error_stack = [] # list of errors on each parameter
best_model_stack = [] # list of best fits to the dats
result_stack = [] # lmfit results
dely_stack = [] # lmfit errors of the fits

for i in range(len(standard_table)):
    # Define initial parameters for Gaussian fitting
    params = Parameters()
    params.add('amp1', value=-1e-5, min=-1, max=0)
    params.add('c1', value=lines_table[sc1_mask][19]['Observed'], min=lines_table[sc1_mask][19]['Observed']-0.5e-4,max=lines_table[sc1_mask][19]['Observed']+0.5e-4)
    params.add('std1', value=spec_res)

    params.add('amp2', value=-1e-5, min=-1, max=0)
    params.add('c2', value=line_center, min=line_center-0.5e-4,max=line_center+0.5e-4)
    params.add('std2', value=spec_res)

    params.add('amp3', value=-1e-5, min=-1, max=0)
    params.add('c3', value=lines_table[si1_mask][2]['Observed'],min=lines_table[si1_mask][2]['Observed']-0.5e-4,max=lines_table[si1_mask][2]['Observed']+0.5e-4)
    params.add('std3', value=spec_res)

    params.add('amp4', value=-1e-5, min=-2e-5, max=0)
    params.add('c4', value=lines_table[sc1_mask][20]['Observed'],min=lines_table[sc1_mask][20]['Observed']-0.5e-4,max=lines_table[sc1_mask][20]['Observed']+0.5e-4)
    params.add('std4', value=spec_res)

    params.add('amp5', value=-1e-5, min=-1, max=0)
    params.add('c5', value=lines_table[na1_mask][1]['Observed'],min=lines_table[na1_mask][1]['Observed']-0.5e-4,max=lines_table[na1_mask][1]['Observed']+0.5e-4)
    params.add('std5', value=spec_res)


    # mathematical/logic constraints on parameters
    # enforce same widths for the Scandium and Sodium lines
    # params['std4'].expr = 'std1'
    # params['std5'].expr = 'std2'

    # Adding a parameter for the wavelength differences
    # params.add('delta', value=1e-3)

    # Constraining the center parameters to have the same difference
    # params['c2'].expr = 'c1 + delta'
    # params['c3'].expr = 'c1 + 2*delta'
    # params['c4'].expr = 'c1 + 3*delta'
    # params['c5'].expr = 'c1 + 4*delta'

    # Define the model
    model = Model(ig.five_gaussian, nan_policy='omit')

    # Perform the fit
    result = model.fit(norm_flux_stack[i], params,
                       x=wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights=1 / raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])

    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit 
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
                         result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
                         result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value,
                         result.params['amp5'].value, result.params['c5'].value, result.params['std5'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
                               result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
                               result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr,
                               result.params['amp5'].stderr, result.params['c5'].stderr, result.params['std5'].stderr])

    best_model_stack.append(result.best_fit)
    # print(result.fit_report())

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = np.array(params_error_stack)[0::3]

centers = params_arr[1::3]
center_err = np.array(params_error_stack)[1::3]

sigmas = params_arr[2::3]
sigmas_err = np.array(params_error_stack)[2::3]

# Generate the fits for each Gaussian component
na1_2_2062_fit = []
sc1_2_2058_fit = []
si1_2_2068_fit = []
sc1_2_2071_fit = []
na1_2_2089_fit = []
num_gauss = 5

for i in range(len(standard_table)):
    sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
    na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
    si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
    sc1_2_2071_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))
    na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[4::num_gauss][i], centers[4::num_gauss][i], sigmas[4::num_gauss][i]))

for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], norm_flux_stack[i], c='black', alpha=0.2, label='Data')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], c='red', lw=1, label='Best Model')
    # ax1.errorbar(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], yerr=dely_stack[i],fmt='ro', ms=3, lw=1, label='Best Model')
    
    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
        best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
        alpha=0.2, color='red')

    # plot the locations of fitted centers
    ax1.scatter(y=[1.03] * len(centers[i * num_gauss:(i + 1) * num_gauss]), x=centers[i * num_gauss:(i + 1) * num_gauss], marker="|", s=100, color='black')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2058_fit[i], ls='-.',
             label=r'$\bf{Sc~I~2.2058}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2062_fit[i], ls='-.',
             label=r'$\bf{Na~I~2.2062}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], si1_2_2068_fit[i], ls='-.',
             label=r'$\bf{Si~I~2.2068}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
            
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2071_fit[i], ls='-.',
             label=r'$\bf{Sc~I~2.2071}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2089_fit[i], ls='-.',
             label=r'$\bf{Na~I~2.2089}$' "\n" rf"A={amps[4::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[4::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[4::num_gauss][i]:.3e}")
    # ax1.set_ylim(top=1.05)
    ax1.set_ylabel('Flux + Constant')

    # Calculate residuals
    residuals = norm_flux_stack[i] - best_model_stack[i]
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i] - np.mean(norm_flux_stack[i]))**2)
    R2 = 1 - (sqsum_res / sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.4f}', fontsize=12)

    plt.suptitle(rf"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
    ax1.legend(bbox_to_anchor=(1, 1), fontsize=12)
    plt.show()
    # print(result_stack[i].fit_report())
In [ ]:
# params_stack = []
# params_error_stack = []

# pcov_stack = []
# best_model_stack = []

# for i in range(len(standard_table)):
#     # Define initial parameters for Gaussian fitting
#     init_params = (-0.05, lines_table[sc1_mask][19]['Observed'], spec_res,
#                           -0.5, line_center, spec_res,
#                           -0.05, lines_table[si1_mask][2]['Observed'], spec_res,
#                           -0.5, lines_table[na1_mask][1]['Observed'], spec_res)

#     # bounds for each parameter
#     # gauss_bounds = ([-1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf],
#     #                 [0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf])

#     # Perform Gaussian fitting for the current source
#     # add error keyword then propogate for equiv width
#     # compare equiv widths between objects of the same spectral types
#     popt, pcov, param_err, best_model = ig.model_fit(ig.four_gaussian,
#                                         wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
#                                         norm_flux_stack[i],
#                                         raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
#                                         init_params,
#                                         maxfev=5000)

#     # Append the best_model, error and parameters stacks
#     params_error_stack.extend(param_err)

#     params_stack.extend(popt)
#     params_arr = np.array(params_stack)

#     # best fit amplitudes
#     amps = params_arr[0::3]
#     amps_err = params_error_stack[0::3]

#     # best fit centers(means)
#     centers = params_arr[1::3]
#     center_err = params_error_stack[1::3]

#     # best fit widths
#     sigmas = params_arr[2::3]
#     sigmas_err = params_error_stack[2::3]

#     pcov_stack.append(pcov)
#     best_model_stack.append(best_model)

# na1_2_2062_fit = []
# sc1_2_2058_fit = []
# si1_2_2068_fit = []
# na1_2_2089_fit = []

# for i in range(len(standard_table)):
#     sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::4][i],centers[0::4][i],sigmas[0::4][i]))
#     na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::4][i],centers[1::4][i],sigmas[1::4][i]))
#     si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[2::4][i],centers[2::4][i],sigmas[2::4][i]))
#     na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[3::4][i],centers[3::4][i],sigmas[3::4][i]))

# for i in range(len(standard_table)):
#     # Create subplots with adjusted size ratios
#     fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black',alpha=0.2, label='Data')
    
#     ax1.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], color='red', label='Best Model')
#     ax1.scatter(y=[1.03] * len(centers[i * 4:(i + 1) * 4]), x=centers[i * 4:(i + 1) * 4], marker="|", s=50, color='black')

#     # for j in range(len(amps[::3])):
#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], sc1_2_2058_fit[i],ls='-.',
#         label= r'Sc I 2.2058 Fit' "\n" rf"A={amps[0::4][i]:.3e}" "\n" rf"$\mu$={centers[0::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2062_fit[i],ls='-.',
#         label= r'Na I 2.2062 Fit' "\n" rf"A={amps[1::4][i]:.3e}" "\n" rf"$\mu$={centers[1::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], si1_2_2068_fit[i],ls='-.',
#         label = r'Si I 2.2068 Fit' "\n" rf"A={amps[2::4][i]:.3e}" "\n" rf"$\mu$={centers[2::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2089_fit[i],ls='-.',
#         label= r'Na I 2.2089 Fit' "\n" rf"A={amps[3::4][i]:.3e}" "\n" rf"$\mu$={centers[3::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::4][i]:.3e}")

#     ax1.set_ylim(top=1.05)
#     ax1.set_ylabel('Flux + Constant')
    
#     # Calculate residuals and R^2
#     residuals = (norm_flux_stack[i] - best_model_stack[i])
#     sqsum_res = np.sum(residuals**2)
#     sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
#     R2 = 1-(sqsum_res/sqsum)

#     # Plot residuals
#     ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
#     ax2.axhline(y=0, color='grey', linestyle='--')
#     # ax2.set_ylim(-0.1,0.1)
#     ax2.set_xlabel(r'Wavelength [$\mu$m]')
#     ax2.set_ylabel('res')
#     ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

#     plt.suptitle(rf"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

#     ax1.legend(bbox_to_anchor=(1,1), fontsize=12)
#     plt.show()
In [ ]:
ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []
ew4_stack = []
ew5_stack = []

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):

    # ew1 = (ig.gaussian_area(amps[0::5][i],sigmas[0::5][i])) # integrate the gaussian
    # ew1_stack.append(ew1)

    # ew2 = (ig.gaussian_area(amps[1::5][i],sigmas[1::5][i]))
    # ew2_stack.append(ew2)

    # ew3 = (ig.gaussian_area(amps[2::5][i],sigmas[2::5][i]))
    # ew3_stack.append(ew3)

    # ew4 = (ig.gaussian_area(amps[3::5][i],sigmas[3::5][i]))
    # ew4_stack.append(ew4)

    # ew5 = (ig.gaussian_area(amps[4::5][i],sigmas[4::5][i]))
    # ew5_stack.append(ew5)
    ew1 = np.trapz(1-na1_2_2062_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
    ew1_stack.append(ew1)

    ew2 = np.trapz(1-sc1_2_2058_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2_stack.append(ew2)

    ew3 = np.trapz(1-si1_2_2068_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3_stack.append(ew3)

    ew4 = np.trapz(1-sc1_2_2071_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew4_stack.append(ew4)

    ew5 = np.trapz(1-na1_2_2089_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew5_stack.append(ew5)

standard_table.loc[:,'ew_na1_2_2062'] = ew1_stack
standard_table.loc[:,'ew_sc1_2_2058'] = ew2_stack
standard_table.loc[:,'ew_si1_2_2068'] = ew3_stack
standard_table.loc[:,'ew_sc1_2_2071'] = ew4_stack
standard_table.loc[:,'ew_na1_2_2089'] = ew5_stack

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2058'],ls='-',marker='o', label=f"Sc I 2.2058 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2062'],ls='-',marker='o', label=f"Na I 2.2062 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_si1_2_2068'],ls='-',marker='o', label=f"Si I 2.2068 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2071'],ls='-',marker='o', label=f"Sc I 2.2071 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2089'],ls='-',marker='o', label=f"Na I 2.2089 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(1e-6,1e-3)
plt.xticks(rotation=45)
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
# plt.yscale('log')

plt.legend(bbox_to_anchor=(1,1))
plt.show()

Mg/Al Interval ~ 2.11¶

Mg 2.1065 = lines_table[mg1_mask][0]

Al 2.1098 = lines_table[al1_mask][0]

In [ ]:
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-160,30),(-80,20),(100,20),(200,20),(270,20),(430,20)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
            plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
    plt.show()
In [ ]:
# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []

# Define initial parameters for Gaussian fitting
init_params = (-0.5, line_center, spec_res,
               -0.5, lines_table[al1_mask]['Observed'][0], spec_res)

for i in range(len(standard_table)):    
    # Define initial parameters for Gaussian fitting
    params = Parameters()

    params.add('amp1', value = 1e-5, max = 0)
    params.add('c1', value = line_center)
    params.add('std1', value = spec_res)

    params.add('amp2', value = 1e-5, max = 0)
    params.add('c2', value = lines_table[al1_mask]['Observed'][0])
    params.add('std2', value = spec_res)

    # params['std1'].expr = 'std2' # constrain widths to be the same

    # Define Model
    model = Model(ig.two_gaussian, nan_policy='omit')

    result = model.fit(norm_flux_stack[i], params,
                       x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr])

    best_model_stack.append(result.best_fit)

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = params_error_stack[0::3]

centers = params_arr[1::3]
center_err = params_error_stack[1::3]

sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]

# Gaussian model fits for each source
mg1_2_1065_fits = []
al1_2_1098_fits = []
num_gauss = 2

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    mg1_2_1065_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))
    al1_2_1098_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))

ew1_stack = []
ew2_stack = []
# ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
    ew1 = np.trapz(1-mg1_2_1065_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-al1_2_1098_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)

standard_table['ew_mg1_2_1065'] = ew1_stack
standard_table['ew_al1_2_1098'] = ew2_stack

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label='Gaussian Fit')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                     best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
                     alpha=0.2, color='red')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_1065_fits[i], ls=':',
     label=r'Mg I 2.1065 $\mu$m' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], al1_2_1098_fits[i], ls=':',
     label=r'Al I 2.1098 $\mu$m' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')

    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1))

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

    plt.suptitle(rf"Mg/Al Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
    plt.show()
In [ ]:
fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_1065'],ls='--',marker='p', label=f"Mg I 2.1065 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_al1_2_1098'],ls='--',marker='p', label=f"Al I 2.1098 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2838'],ls='--',marker='p', label=f"Fe I 2.2838$\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(-0.0001,0.0002)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1,1))
plt.show()

Ti Region (+ Fe)¶

In [ ]:
# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][85]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-600,30),(-350,50),(-150,60),(130,10),(400,100)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
params_stack = []
params_error_stack = []

pcov_stack = []
best_model_stack = []

# Define initial parameters for Gaussian fitting
init_params = (-0.2, lines_table[ti1_mask]['Observed'][44], spec_res,
               -0.2, lines_table[ti1_mask]['Observed'][45], spec_res,
                -0.2, line_center, spec_res,
                -0.1, lines_table[fe1_mask]['Observed'][86], spec_res,
                -0.1, lines_table[ti1_mask]['Observed'][47], spec_res)

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.five_gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i],
                                        raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)
    # Append the best_model, error and parameters stacks
    params_error_stack.extend(param_err)

    params_stack.extend(popt)
    params_arr = np.array(params_stack)

    amps = params_arr[0::3]
    amps_err = params_error_stack[0::3]

    centers = params_arr[1::3]
    center_err = params_error_stack[1::3]

    sigmas = params_arr[2::3]
    sigmas_err = params_error_stack[2::3]

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

ti1_2_2217_fit = []
ti1_2_2239_fit = []
fe1_2_2263_fit = [] # maybe fix wavelength offsets
fe1_2_2266_fit = []
ti1_2_2280_fit = []
num_gauss = 5

for i in range(len(standard_table)):
    ti1_2_2217_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))

    ti1_2_2239_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))

    fe1_2_2263_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[2::num_gauss][i],centers[2::num_gauss][i],sigmas[2::num_gauss][i]))

    fe1_2_2266_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[3::num_gauss][i],centers[3::num_gauss][i],sigmas[3::num_gauss][i]))

    ti1_2_2280_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[4::num_gauss][i],centers[4::num_gauss][i],sigmas[4::num_gauss][i]))
    

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], alpha=0.2, c='black', label='Data')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', lw=2, label='Best Model')
    
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2217_fit[i], ls='--',
        label = 'Ti I 2.2217 Fit' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2239_fit[i], ls='--',
        label = 'Ti I 2.2239 Fit' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2263_fit[i], ls='--',
        label = 'Fe I 2.2263 Fit' '\n' rf'A={amps[2::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[2::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[2::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2266_fit[i], ls='--',
        label = 'Fe I 2.2266 Fit' '\n' rf'A={amps[3::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[3::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[3::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2280_fit[i], ls='--',
        label = 'Ti I 2.2280 Fit' '\n' rf'A={amps[4::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[4::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[4::num_gauss][i]:.3e}')

 
    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1), fontsize=12)

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

    plt.suptitle(rf"Ti Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    ew1 = np.trapz(1-fe1_2_2263_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
    ew1_stack.append(ew1)

    ew2 = np.trapz(1-fe1_2_2266_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2_stack.append(ew2)

    ew3 = np.trapz(1-ti1_2_2280_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3_stack.append(ew3)

standard_table.loc[:,'ew_fe1_2_2263'] = ew1_stack
standard_table.loc[:,'ew_fe1_2_2266'] = ew2_stack

standard_table.loc[:,'ew_ti1_2_2280'] = ew3_stack

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-',marker='o', label=f"Fe I 2.2263 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2266']),ls='-',marker='o', label=f"Fe I 2.2266 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_ti1_2_2280']),ls='-',marker='o', label=f"Fe I 2.2280 $\mu$m")

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")


plt.legend(bbox_to_anchor=(1, 1))
plt.show()

Ca Region¶

Ca I 2.2614

Fe I 2.2626 = lines_table[fe1_mask][29]

Ca I 2.2631

Ca I 2.2657

In [ ]:
# import matplotlib
# matplotlib.use('TkAgg')  # Set the backend to TkAgg for interactivity
# # Define the region for fitting
# line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
# line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength

# # from igrins_mod import local_continuum_fit
# continuum_stack = []
# norm_flux_stack = []
# reg_idx_stack = []

# # Helper function to select regions interactively
# def select_regions(wavelengths, flux, num_regions):
#     plt.plot(wavelengths, flux)
#     plt.title('Click to select points for region exclusion')
#     plt.xlabel('Wavelength')
#     plt.ylabel('Flux')

#     # Select points interactively
#     selected_points = plt.ginput(n=num_regions*2, timeout=0) # n=num_regions*2 because we need start and end points for each region
#     plt.close()

#     regions = []
#     for i in range(0, len(selected_points), 2):
#         start = selected_points[i][0]
#         end = selected_points[i+1][0]
#         regions.append((start, end))
    
#     return regions

# # number of regions to select interactively
# num_regions = 4
# poly_deg = 3

# for i in range(len(standard_table)):
#     idx1 = np.searchsorted(wavelen_stack[:, i], line_center ) - 150
#     idx2 = np.searchsorted(wavelen_stack[:, i], lines_table[ca1_mask]['Observed'][2]) + 150

#     wavelengths = wavelen_stack[idx1:idx2, i]
#     flux = raw_flux_stack[idx1:idx2, i]
    
#     # Select regions interactively
#     regions = select_regions(wavelengths, flux, num_regions)
    
#     # Prepare the regions for the local_continuum_fit function
#     regions = [(int(np.searchsorted(wavelengths, start)), int(np.searchsorted(wavelengths, end))) for start, end in regions]
    
#     fig = plt.figure(figsize=(15, 3))
#     continuum, region_indices = ig.local_continuum_fit(wavelengths,
#                                                        flux,
#                                                        poly_order=poly_deg,
#                                                        line_center=line_center,
#                                                        spec_res=spec_res,
#                                                        regions=regions)
#     # Append indices to list of indices
#     reg_idx_stack.append(region_indices)
    
#     # Append to list of the local continuum arrays
#     continuum_stack.append(continuum)

#     # Normalize flux by dividing the flux by the continuum and create a list
#     norm_flux = flux[region_indices[0][0]:region_indices[-1][1]] / continuum_stack[i]
#     norm_flux_stack.append(norm_flux)

#     plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], flux[region_indices[0][0]:region_indices[-1][1]])
#     plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], continuum_stack[i])

#     for j in range(len(region_indices)):
#         plt.axvspan(wavelengths[region_indices[j][0]], wavelengths[region_indices[j][1]], color='red', alpha=.2)

#     plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

#     plt.show()
In [ ]:
# Define the region for fitting
line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-100,50), (50,30), (240,50),(550,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
np.abs(lines_table[ca1_mask]['Observed'][0] - lines_table[ca1_mask]['Observed'][1])
# lines_table[ca1_mask]['Observed'][1]
# lines_table[ca1_mask]['Observed'][2]
# lines_table[fe1_mask]['Observed'][104]
Out[ ]:
0.0017000000000000348
In [ ]:
# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []

# Define initial parameters for Gaussian fitting
# init_params = (-0.5, lines_table[ca1_mask]['Observed'][0], spec_res,
#                -0.5, lines_table[ca1_mask]['Observed'][1], spec_res,
#                -0.5, lines_table[ca1_mask]['Observed'][2], spec_res,
#                -0.5, lines_table[fe1_mask]['Observed'][104], spec_res)
from scipy.signal import savgol_filter



for i in range(len(standard_table)):    
    # Define initial parameters for Gaussian fitting
    params = Parameters()

    params.add('amp1', value = 1e-5, max = 0)
    params.add('c1', value = lines_table[ca1_mask]['Observed'][0])
    params.add('std1', value = spec_res)

    params.add('amp2', value = 1e-5, max = 0)
    params.add('c2', value = lines_table[ca1_mask]['Observed'][1])
    params.add('std2', value = spec_res)

    params.add('amp3', value = 1e-5, max = 0)
    params.add('c3', value = lines_table[ca1_mask]['Observed'][2])
    params.add('std3', value = spec_res)

    params.add('amp4', value = 1e-5, max = 0)
    params.add('c4', value = lines_table[fe1_mask]['Observed'][104])
    params.add('std4', value = spec_res)

    params['std1'].expr = 'std2' # constrain widths to be the same
    params['std1'].expr = 'std3'
    # params['std1'].expr = 'std2'

    # Define Model
    model = Model(ig.four_gaussian, nan_policy='omit')

    result = model.fit(savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), params,
                       x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])

    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
                         result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
                         result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
                               result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
                               result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])

    best_model_stack.append(result.best_fit)

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = params_error_stack[0::3]

centers = params_arr[1::3]
center_err = params_error_stack[1::3]

sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]

# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
    ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
    ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
    fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 5 \sigma$')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), c='black', alpha=0.5)

    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                    best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
                    alpha=0.2, color = 'red' )

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
            label=r'Ca I 2.2614 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
        label=r'Ca I 2.2631 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.',
        label=r'Ca I 2.2657 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.', c='magenta',
        label=r'Fe I 2.2626 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    
    # plt.axvline()
    plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1))

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
    
    plt.show()
In [ ]:
ew1_stack = []
ew2_stack = []
ew3_stack = []
ew4_stack = []

for i in range(len(standard_table)):
    ew1 = np.trapz(1-ca1_2_2614_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-ca1_2_2631_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3 = np.trapz(1-ca1_2_2657_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew4 = np.trapz(1-fe1_2_2626_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)
    ew3_stack.append(ew3)
    ew4_stack.append(ew4)

standard_table['ew_ca1_2_2614'] = ew1_stack
standard_table['ew_ca1_2_2631'] = ew2_stack
standard_table['ew_ca1_2_2657'] = ew3_stack
standard_table['ew_fe1_2_2626'] = ew4_stack

fig = plt.figure(figsize=(15,5))

# plt.plot(standard_table['Spectral_Type'], standard_table['ew_sc1_2_2058'],ls='--',marker='p', label=f"Sc I 2.2058 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2062'],ls='--',marker='p', label=f"Na I 2.2062 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_si1_2_2068'],ls='--',marker='p', label=f"Si I 2.2068 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2089'],ls='--',marker='p', label=f"Na I 2.2089 $\mu$m")

# plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-.',marker='o', label=f"Fe I 2.2263 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2266'],ls='-.',marker='o', label=f"Fe I 2.2266 $\mu$m")

plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2614'],ls=':',marker='s', label=f"Ca I 2.2614 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2631'],ls=':',marker='s', label=f"Ca I 2.2631 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2657'],ls=':',marker='s', label=f"Ca I 2.2657 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2626'],ls=':',marker='s', label=f"Fe I 2.2626 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(-0.0001,0.0005)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1, 1))
plt.show()

Mg, Ca, Fe ~ 2.28 micron¶

Mg I 2.2813 = lines_table[mg1_mask][8]

Fe I 2.2838 = lines_table[fe1_mask][115] (possible 116)

Ca ? = 2.285?

In [ ]:
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][8]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-150,10),(-100,10),(70,10),(300,10),(500,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
In [ ]:
# min_wave = line_center - (100 * spec_res)
# max_wave = line_center + (100 * spec_res)

# Optimal parameters that will be appended from scipy curve fit below
amplitude1_stack = []
center1_stack = []
sigma1_stack = []

amplitude2_stack = []
centeR2_stack = []
sigma2_stack = []

amplitude3_stack = []
center3_stack = []
sigma3_stack = []

pcov_stack = [] # Covariant matrices
best_model_stack = [] # Best model fits

flux_constant = np.linspace(0,-1,len(standard_list))

# Define initial parameters for Gaussian fitting
init_params = ( -0.5, line_center, spec_res,
                -0.5, lines_table[fe1_mask]['Observed'][114], spec_res,
                -0.5, 2.2825, spec_res)

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.three_gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)
    # Append the Gaussian parameters, covariant matrices, and model fits as we loop through
    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma1_stack.append(popt[2])

    amplitude2_stack.append(popt[3])
    centeR2_stack.append(popt[4])
    sigma2_stack.append(popt[5])

    amplitude3_stack.append(popt[6])
    center3_stack.append(popt[7])
    sigma3_stack.append(popt[8])

    # amplitude4_stack.append(popt[7])
    # center4_stack.append(popt[8])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# Gaussian model fits for each source
mg1_2_2813_fits = []
fe1_2_2818_fits = []
ca1_2_2825_fits = []

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    mg1_2_2813_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma1_stack[i]))
    fe1_2_2818_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude2_stack[i],centeR2_stack[i],sigma2_stack[i]))
    ca1_2_2825_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude3_stack[i],center3_stack[i],sigma3_stack[i]))

ew1_stack = []
ew2_stack = []
ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
    ew1 = np.trapz(1-mg1_2_2813_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-fe1_2_2818_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3 = np.trapz(1-ca1_2_2825_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)
    ew3_stack.append(ew3)


standard_table['ew_mg1_2_2813'] = ew1_stack
standard_table['ew_fe1_2_2818'] = ew2_stack
standard_table['ew_ca1_2_2825'] = ew3_stack

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red',label='Gaussian Fit')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_2813_fits[i],ls=':',label=r'Mg I 2.2813 $\mu$m')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2818_fits[i],ls=':',label=r'Fe I 2.2818 $\mu$m')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m')
    
    # plt.axvline()
    ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    ax1.set_ylabel('Flux + Constant')
    ax1.set_xlabel(r'Wavelength [$\mu$m]')
    ax1.legend(bbox_to_anchor=(1,1))
    
    plt.show()
c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[40], line 95
     92 ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m')
     94 # plt.axvline()
---> 95 ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
     97 ax1.set_ylabel('Flux + Constant')
     98 ax1.set_xlabel(r'Wavelength [$\mu$m]')

AttributeError: 'Axes' object has no attribute 'suptitle'
In [ ]:
fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_2813'],ls='--',marker='p', label=f"Mg I 2.2813 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2818'],ls='--',marker='p', label=f"Fe I 2.2818 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2825'],ls='--',marker='p', label=f"Ca I 2.2825 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

plt.ylim(-0.0001,0.0002)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1,1))
plt.show()
In [ ]:
# standard_table.to_csv('standard_table_v2.txt')